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Abstract 



We use our recently proposed accelerated dynamics algorithm (Tiwary and 



van de Walle, 2011) to calculate temperature and stress dependence of ac- 
tivation free energy for surface nucleation of dislocations in pristine Gold 
nanopillars under realistic loads. While maintaining fully atomistic resolu- 
tion, we achieve the fraction of a second time-scale regime. We find that the 
activation free energy depends significantly and non-linearly on the driving 
force (stress or strain) and temperature, leading to very high activation en- 
tropies. We also perform compression tests on Gold nanopillars for strain 
rates varying between 7 orders of magnitudes, reaching as low as 10^/s. Our 
calculations bring out the perils of high strain-rate Molecular Dynamics cal- 
culations: we find that while the failure mechanism for (001) compression 
of Gold nanopillars remains the same across the entire strain-rate range, the 
elastic limit (defined as stress for nucleation of the first dislocation) depends 
significantly on the strain-rate. We also propose a new methodology that 
overcomes some of the limits in our original accelerated dynamics scheme 
(and accelerated dynamics methods in general). We lay out our methods in 
sufficient details so as to be used for understanding and predicting deforma- 
tion mechanism under realistic driving forces for various problems. 
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1. Introduction 



Forming a correct picture of dislocation nucleation is central to our under- 
standing of deformation mechanisms at the nano-scale. The initial discoveries 



by Uchic and subsequent work by various groups ( 


Uchic et al. 


2004t Wu et al. 


2005 


Greer et al. 


2005 


20091 Volkert and Lilleodden 


2006[ S 


lan et al. 


2007 


Jennings et al. 


2010 


I 


A et aL 


2010| Jennings et al. , 


2011 


) have now estab- 



lished that there is a marked increase of yield strength as the specimen size 
decreases, with significant strain-rate dependence as well. These observations 
have generally been attributed to the scarcity of dislocation sources (such as 
Frank-Read sources) in nanosized samples, and having to nucleate disloca- 



tions in a perfect crystal (perfect apart from presence of surfaces) ( Greer and 



Nix 


2006 


Nix et al. 


2007 ZhuandLi 


2010 



such there have been numerous attempts to link simulations of dislocation nu- 
cleation processes to experimentally observed mechanical behavior - in fact a 



lot of crucial insight has come from simulations (Bulatov et al. 



et al. 



et al. 



1997 Moriarty et al. 



2007 Zhu et al. 



(Landman et al. 1990 



2008 



1997 



2002| |Li et aL] [2002| [Abraham et al.] [2002| [Cao 



Marian 



Ryu et al. 2011 ). Nanoindentation experiments 
Schuh et al. , 2005 Schuh, 2006), Scanning Electron 
Microscopy combined with Nanoindentation (Jang and Greer 2010), and 



High Resolution Transmission Electron Microscopy (HRTEM) ( Zheng et al 



2010) are now sufficiently advanced for one to hope for a direct match be- 



tween simulations and experiments ( 


Van Vliet et al. 


2003 


Zhao et al. 


2003 


Greer et al. , 


2008 


Zheng et al. 


2010 


)• 



There are three broad classes of techniques for such simulations, often 
used in conjunction with each other and along with approaches such as Tran- 
sition State Theory: classical Molecular Dynamics, continuum based meth- 
ods and ab initio techniques. Ab initio simulations, though they often provide 



insight into mechanical behavior 


Arias and Joannopoulos 


1994 




Sitch et al. 


1995 


Woodward and Rao, 


2002 


Woodward et al. 


2008 


Carter 


2008), are 



still restricted to very small sizes, less than a hundred atoms typically (al- 
though there are several promising attempts at bridging this length-scale gap 



Fago et al. , 


2004 


Gavini et al. 


2007a|b 


Trinkle 



2008)). The achievable time-scales are also typically restricted to less than 
a few picoseconds. Continuum based methods are another elegant option 
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capable of dealing with a variety of length and time scales, though they suf- 
fer from not providing atomic scale resolution and assuming elastic behavior 



even at dislocation cores (Hirthe and Lothe 1982 Cai et al. 2006 Wein 



berger et al. 2012). Classical Molecular Dynamics (MD) can be helpful in 



gaining quantitative insight into mechanical behavior at various length scales 



(nanometers to microns or larger) ( Zhou et al. , 1998 Chang et al 



et al. 2006 Park et al. 2006 Rabkin et al. 2007 Li et al. 2010). MD does 



2002; Diao 



not assume much apart from the form of interatomic interaction, which is 
typically developed by fitting to first principles or experimental data. The 



availability of quality interatomic force- fields (Foiles et al. , 1986 Justo et al. 



1998 Grochola et al. 2005; Rabkin et al. , 2007) and increase in computer 



power has led to a tremendous increase in the popularity of MD over the last 
decade. 



However, most of the interesting dynamics happens only as the system 
moves from one energy basin to another through infrequent, rare events. 
Most of the simulation time gets spent with the system staying stuck in 
some energy basin (Laio and Parrinello, 2002). This behavior, combined 



with the femtosecond timestep required for total energy conservation, gives 



rise to a major limitation of MD: the time-scale problem (Voter et al. 2002 



Kushima et al. 2011 Li et al. 2011 ). Even with the advent of powerful super- 



computers, MD simulations are unable to reach more than a few nanoseconds 
of time if the system size is more than a few hundred atoms. Thus while lab- 
oratory strain-rates are typically in the range 10~^-10'^/sec, with correspond- 
ing activation free energies being around SO/c^T, MD is unable to go slower 
than 10^/sec strain-rate, corresponding to free energies of around Sfc^T or 



Walle 



lower(|Zhu et al. 2008 Hara and Li 2010 Ryu et al. 2011 Tiwary and van de 
2011[ ). One approach to get around this shortcoming is to perform 



temperature Nudged Elastic Band calculation (Henkelman et al. 2000) of 



the activation free energy and how it varies with applied stress, and then 
either assume it to be temperature independent, or assume a phenomenolog- 
ical model for its variation with temperature (such as multiplying it with an 



empirical temperature dependent scaling factor) (Weinberger et al. , 2012) 



These approaches can sometimes work well, but as shown by Hara and Li 



(2010); Ryu et al. (2011) and in this current work, can sometimes lead to 



significant inaccuracies in the predictions, such as errors of several orders of 
magnitude in the nucleation rate, or even qualitatively incorrect phase tran- 



sitions (Warner and Curtin, 2009). 
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We recently proposed a hybrid method that combines the strengths of 
MD and Monte Carlo (MC) simulations in an easily implementable manner 
and enables one to reach milliseconds and longer times for several thousands 



of atoms (Tiwary and van de Walle 2011 ). The method was recently applied 



to vacancy diffusion in BCC Fe at low temperature as well as calculation 
of stress-strain behavior for Au nanowire compression. In both cases it was 
found to predict correct dynamics, with excellent scaling in computational 
efficiency with system size. The algorithm was particularly designed to be 
used on massively parallel computer systems. The main advantages of this 
new method over existing accelerated MD techniques (see, for e.g.. Voter 



( [1M7| ); [Voter et al.| ( [2002| ); |Miron and Fichthorn] ( |2003| ) ; |Lu et al.| ( [2010| )) are 
that (i) it provides a statistically more accurate "real" time scale (which is 
important when determining the actual strain rate in a simulation and time- 
dependent forces in general), (ii) it does not rely on transition state theory 
(which is important when the object of interest is the entropy of activation or 
migration), (iii) it does not require the specification of the degrees of freedom 
of interest (which is crucial when the mechanisms are complex and involve 
the movement of many atoms) and (iv) the efficiency of estimating the "real" 
time scale (as in (i)) improves linearly with number of computer processors 
employed for the calculation. 

It has remained an unsolved problem so far to design and peform fully 
atomistic simulations that could provide a picture of temperature depen- 
dent activation free energies of dislocation nucleation from surfaces at real- 
istic loads and loading rates. Such a picture is key to linking experimental 



results with simulation predictions (Warner and Curtin, 2009 Ryu et al 



2011 Weinberger et al. , 2012). The critical nucleus for surface nucleation 
can be as small as a few atomic planes, thus questioning the applicabil- 
ity of continuum methods. As for classical MD, the time-scale achievable 
is several orders of magnitude smaller than experiments, thus limiting MD 
simulations to regimes of extremely high nucleation rates. With our recently 
proposed hybrid MC-MD method that allows us to achieve extended time- 
scales while still maintaining atomistic resolution, we are able to study the 
temperature dependence of activation parameters for surface nucleation of 
dislocations in pristine nanowires and obtain several significant results in an 
activation regime actually achievable in laboratory experiments. The spe- 
cific problem we consider pertains to several nano-indentation experiments 
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where it was found that even if the apphed stress on a sample is in the elastic 
regime, yielding could occur after a certain statistically distributed waiting 



time (Zuo et al. 2005 Zuo and Ngan, 2006 Mason et al. 2006). We per- 



form fully atomistic simulations of this time-dependent incipient plasticity 
behavior in Gold nanowires, reaching hundreds of milliseconds time-scales for 
several thousand atoms. After collecting statistics for various temperatures 
and applied stresses, we then derive the full picture of stress and temperature 
dependence of the activation free energy. 

In Section 2, we provide a brief summary of out hybrid method in or- 
der to facilitate further discussion and keep this paper self-contained. For a 



more detailed explanation of terms we refer the reader to Tiwary and van de 



Walle (2011). In the current paper, we also propose a new adiabatic switch- 



ing technique that significantly reduces the number of input parameters in 
our hybrid MC-MD approach and eliminates some of the fundamental limi- 
tations of our earlier implementation (that were shared by related algorithms 



(Voter, 1997)). The algorithms employed here make it possible to achieve 
linear scaling in efficiency of estimating the accelerated time as the number 
of parallel processors employed is increased. We describe our algorithm and 
its implementation in sufficient detail for researchers to be able to use it for 
their problem of interest, and hope that it will be found helpful for modeling 
a variety of mechanical behavior problems. 



2. Details of calculations 



2.1. Choice of interatomic potential 

There are several good potentials available for modeling mechanical be- 



havior of Gold (Foiles et al. 1986 Cai and Ye, 1996; Park and Zimmerman 



2005 Grochola et al. , 2005). The embedded atom method potential by Gro- 



chola et al. (2005) gives very realistic values for the surface energy and the 



stacking fault energy (Rabkin et al. , 2007): the stacking fault energy from 
the potential by Grochola et al. (2005) is 42 mJ/m^ while the experimental 



value for it is in the range 32-46 mJ/m^. Since the current paper deals with 
nucleation of dislocations from surfaces, we choose the Grochola potential. 
This potential was also used and found to perform very well in a recently 
published joint computational and experimental work studying dislocation 



behavior in sub-lOnm Gold nanowires (Zheng et al. 2010). Rabkin et al. 
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(2007) provide a critical comparison of this potential with other available 
potentials for various physical properties relevant to the current work. 



2.2. Hybrid stochastic and deterministic technique for achieving realistic time- 
scales 

2.2.1. Summary of ideas 

We recently proposed using a combination of MD and MC techniques for 
achieving long time scales (Tiwary and van de Walle, 2011 ). Our approach is 
built upon minimizing the MD time spent in low-lying energy basins, and in- 
stead using 2 kinds of MC simulations. One (a) seeks to properly thermalize 
the system between infrequent events, thereby minimizing artifical correla- 
tions, and the other (b) provides independent control over the accuracy of 
the time-scale correction. When the potential energy V{x) of the system 
(where a; is a point in the 3-N dimensional configuration space for a system 
with N particles) is above a certain Vq, the system evolves as per regular MD 



see Fig. 1 in Tiwary and van de Walle (2011)). This high energy region of 



the phase space is the one containing the interesting but infrequent events. 
When the system potential energy goes below Vq, we allow MD to continue 
until the system has lost memory of how it entered this well (defined as all 
points X such that V{x) < Vq). We found that a simple and appropriate cri- 
terion to check for this memory loss is when the energy reaches the system's 
mean energy at that temperature (although other choices are possible, such 
as letting MD continue for a sufficiently long, user-specified, time or for a 
random length of time drawn from a user-specified exponential distribution). 
During this thermalization time, the system may escape the well, in which 
case the system simply keeps evolving via MD. Most likely, however, the sys- 
tem will not escape the well during that time. When this happens, we stop 
MD and launch the first MC simulation (called MC a). MC a runs with a 
perfectly uniform potential inside the well, rejecting all moves that lead to 
V{x) > Vq. The purpose of MC simulation a is to generate a new, properly 
thermalized, starting point for MD. MD resumes with positions drawn from 
the last MC a step that visited the boundary of the well(defined rigorously 
in Eq.([2])), and velocities drawn from a Maxwell-Boltzmann distribution in 
the half space pointing outward of the well. Vq can be picked high or low 
depending on the speed-up relative to MD we seek for a particular apphca- 
tion. The method is formally correct for any choice of Vq; a higher choice of 
Vq limits our ability to monitor the detailed dynamics of some events. 
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We also need to estimate the expected value of the time the system would 
have spent in the energy well W, which can be calculated as the reciprocal 
of the flux exiting the well: 



ty^ = lim((- l(a; G S^)))'^ 



(1) 



where the average (■ ■ ■) is taken over x drawn from the well W with a prob- 
ability density proportional to e"^*^^^/'-'^'^^''. ks is Boltzmann's constant, T 
is the temperature and the following definitions hold: 1{A) equals 1 if the 
event A is true and otherwise, is a shell of width w at the boundary of 
the well W, which can be defined in the limit of small w as 



{x : \V{x) - Vol < w\VV{x)\/2} 



(2) 



V denotes the mean projection of a Maxwell-Boltzmann-distributed velocity 
along the unit vector u parallel to W{x), conditional on v ■ u > 0. When 
all atoms have the same mass m,v= ^JksT /2TTm (a general expression can 
be found in [Tiwary and van de Walle (2011)). 



Since Eq.Q involves an average, it can be approximated using MC sim- 
ulations. We make use of the system's ergodicity, replacing the time average 
(that would require us to wait for long times for it to converge) by an en- 
semble average. Thus in parallel to MC a, we launch several instances (as 
many as number of available processors) of a second kind of MC simulation, 
called MC 6, to estimate the time-scale correction. A most straightforward 
implementation of this still won't be as effective in estimating the average in 
Eq.Q because the shell would be visited very rarely. Thus, to improve 
the efficiency in estimating Eq.([T]), we proposed in Tiwary and van de Walle 
( |20ll| ) using a biased potential V*{x), which is the same as the true potential 



V{x) in the high energy regions but lifted up in the energy basins (Voter 



use in this (Voter et al. 



1997 Voter et al.[ 2002|). Several lifting (or biasing) schemes are available for 

as 



2002). A simple importance sampling expression 



detailed in Tiwary anc 
time-correction: 



van de Walle] ( |20ll| )) can thus give us the following 



lim 



-P{V{x)-V*(x))Y 



-^0 (£e-/^(^(-)-^*(-))l(x e S^))* 

where (■■■)* denote expectations taken under a density proportional to 
-pv*{x)^ in which (3 is 1/(A;bT). This approach works well, but there is a 
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fundamental trade-off that limits its usefulness. Lifting the biased potential 
more and more leads to the energy shell being visited more frequently 
and should lead to greater computational efficiency in estimating Eq.Q. 
However a compensating effect leads to a decrease in efficiency beyond a 
certain amount of biasing. This is because increased biasing of the potential 
leads to noisier statistical everaging of the time in Eqs. ([T| and (fsl) - as biasing 
increases, {V{x) — V*{x)) becomes a large number causing e~^^^^~'^*^^''^ to 



dramatically increase. This point has been discussed in detail in 


Miron and 


Fichthorn 


(2003 


); 


Hara and Li 


(2010 


); 


Tiwary and van de Walle 




Mil 


). Note 



that need not be computed for every well the system visits. The overall 
time scaling only the depends on the average escape time, which may be well 
estimated by a converged average of over a number of wells, though not 
necessarily all wells, but instead over a randomly chosen sampling of all the 
wells. 



2.2.2. Adiabatic switching technique 

We now propose a technique that bears some resemblance to adiabatic 



switching methods (see, e.g. (Tuckerman, 2010)) that helps us deal with the 



trade-off discussed above, and also eliminates the need for picking a partic- 
ular biasing scheme. The motivation here is to avoid the statistical noise 
m Eq.(|3|) that arises as the biased potential V*{x) becomes increasingly dif- 
ferent from the true potential V{x). To avoid this noise in sampling, the 
system is continuously, adiahatically switched from V (x) (the true poten- 
tial) to Vo (a fiat potential within the well, identical to the potential used 
in MC simulation a for thermalization). We now formally derive the method. 

Let V (x, a) smoothly interpolate between V (x, 0) = (x) and V (x, 1) = 
Vq. Then we can express the ensemble average in Eq.([T]) as below (working 
in terms of rate = 1/time): 



rate = lim 

ui— s>0 



lim 

ui->0 



/|l(x£5^)e-^^("'°)^x 
/ e-/'^(^'0)c/x 

/e-/3^(^'i)rfx I Je-/3^(^'0)dx 



lim / ^ ^""f g-/3(v-(^,o)-v(^,i)) \ ^ 
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where dx denotes a differential volume in 3-N dimensional configuration space 
for N particles, the integration being performed over entire configuration 
space within the well W and the expected value (■ ■ ■)„ in Eq.(|4]) is defined 
by 

_ /(■■■) e-/^^(---)rfx 

Below we define the term R in Eq. Q and re-express it in a computationally 
tractable form: 



R 




exp (in [ e-^^(^'^)rfx - In / e-^^(^'°)rfx 
exp 

With this we can now write the rate in Eq.Q as 

rate = ^^^v(^-^^'^^ exp (-/^f ( ^^g^^ ) 

" (7) 

We now make a few observations regarding the above expression. It 
involves 3 independent parts. The first is which we already know as 
yJk^Tj2wm for identical atoms. We could keep v inside the ensemble aver- 
age to cover the general case of unequal masses in which v may depend on x. 
The second part in Eq.(7) is lim^_,o (^Mi£§^e-^(^(^'°)-^(^'^))^^. This is non- 

only when x G S"^, and whenever it is non-0, the difference V (x, 0) — V (x, 1) 
is a very small number (see Eq.(|2]))). Since this average is calculated with a 
fiat potential V"(x,l), the boundary x e 5^ is visited frequently, and thus 
the second term in Eq. ([T]) can be evaluated very quickly - in a few MC passes 
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200 400 600 
Number of MC passes 



800 



0.41 



12 3 4 

w (Angstroms) xlO"^ 



Figure 1: (a) The 2nd part in Eq.(7), i.e. ^ U^es^) e-/3(v-(x.o)-v-(x.i)) \^ ^ be eval- 
uated in a very small number of MCT passes as explained in the text, (b) Calculating 
limuj^o \ ^^^^^^^^^e^^^^'^^'"-'^^'^^'^-') ) using linear extrapolation. Calculations are for a 



Au nanowire with 2016 atoms, 2.5nm in diameter and 7.5nm in height, at 300 K. 



as shown in figure la). We calculate M = (^liEi^e-^i^^'''^^-^^^'^^)^ for 
few values of and simple linear extrapolation gives the desired limit, as 
shown in figure 1 ^b). The third part in Eq.(7) is exp (^~f^ Jq ^ ^^1"^'"^ \ da 



da 



Here, the average 



/ dV{x,a) \ 



Qa ~)a does not contain any exponentials, and thus no 
terms that could blow-up and lead to noisy estimates and slow convergence. 



We now need to pick up a switching scheme for V{x, a), i.e. an interpo- 
lation scheme between V{x,0) and V{x, 1). We picked the simplest scheme 
- a linear switching model - and found it to work very well: 

V{x,a) = {1 -a)V{x) + aVo (8) 

With this, Eq.Q for the rate of escaping energy basins bounded by 
V{x) < Vq becomes 



rate = \im v ( ^ ^ g^/^(v-(x,o)-yM)) 



exp j {V{x) - Vo)o.da 



(9) 

Thus to summarize till this point, to calculate Eq.(|9]), we first do a quick 
MC simulation using a fiat potential to get lim^^o ( MiE|^e-^(^(^'°)-^(^'^)) ^ , 
as shown in figure [T] We then vary a adiabatically during the simulation, 
going from a = to a = 1. We perform a series of MC simulations, with 
the Hamiltonian of the system evolving as per Eq.pj) as the simulation time 
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progresses. A typical evaluation of Eq. ^ done as per this scheme is shown 
in figure [2| where we show the change in the following two as a function of 



a: a 



1 / dV{x,a) 



da 



da), related to the 3rd part in Eq.(9), and (b) the 



expected value of the time spent in the energy well W. 




Figure 2: (a) Change in (^—(3 y '^^ia'"'' / '^'^J ^ ^ function of a as the simulation 
progresses (see Eq.([7])). (b) Expected value of the time spent in the energy well W. 
Calculations are for a Au nanowire with 2016 atoms, 2.5nm in diameter and 7.5nm in 
height, at 300 K. 



It may be useful to emphasize one point: The reason we do not immedi- 
ately start MC as soon as V{x) falls below Vq is because the rate expression 
(Equation ([T])) is only valid conditional on the system being initialized at a 
Boltzman-distributed random position within the well. If we start MC at 
the boundary of the well, this assumption is violated. Although it may seem 
that, by running MD within the well for some time before calculating the 
escape rate, we slightly overestimate the time spent in the well, this is not 
the case. The distribution of escape time from the well is independent of the 
time already spent in the well (since it follows an exponential distribution). 
Another way to see that there is no time over-couting is to observe that, dur- 
ing this MD trajectory in the well, there is also a small probability that the 
system escapes the well, so we are not artificially constraining the system to 
remain in the well for a longer time. The above scheme is a new, yet formally 
correct, way to deal with so-called re-crossing events that typically affect the 



accuracy of transition-theory based estimates of escape times (Voter et al. 



2002 Lu et al. 2010). 
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Figure 3: Perspective view of the nanopillar (a) before application of strain and (b) after 
application of strain when partial dislocation has nucleated. Coloring is as per Common 
Neighbor Analysis! Honey cutt and Andersen 1987| ), where green denotes FCC, red denotes 
HCP. In (b), the surface atoms (identified as atoms that are neither FCC nor HCP) have 
been removed to bring the slip plane into clarity. Visualization was carried using the 
package OVITOptuk^skil [20lo|. 



2.3. Simulation set-up and compression testing 

We first report the stress-strain plots for (001) compression of pristine 
cylindrical Au nanowires. The cylinder was initially carved out from perfect 
FCC lattice and before compression, it was 2.5nm in diameter and 7.5nm 
in height, comprising 2016 atoms (see figure |3]) with periodic boundary con- 
ditions imposed along all three directions. The cylinder axis z is also the 
compression axis (001). Thus all sites along the length of the wire are now 
equivalent sites for nucleation. For the other two directions, we do not strictly 
need periodic boundary conditions, but we nevertheless apply it for compu- 
tational ease. The dimensions of the supercell in the x and y directions are 
both around 75A , which is much larger than the range of the EAM potential 
employed (5.5A)(Grochola et al. , 2005). As such there is no artifact from 
the pillar interacting with its images in these two directions. The cylinder 
was first equilibrated for 500 ps before beginning the compression, which 
was carried out by uniformly re-scaling the z-coordinates of all atoms. The 



atomic virial stress was used to obtain the Cauchy stress (Rabkin et al. , 2007) 



The stress at zero nominal strain is non-zero and tensile, and arises from the 
surface stress 



see 



Rabkin et al. (2007); Weinberger et al. (2012) for a more 
detailed explanation). We adjust for this, and as such figure |4| provides the 
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stress span, i.e. the stress at a strain e relative to the stress at strain. We 
present the resulting stress((j) versus nominal strain(e) plots for 2 different 
strain rates e: 5xl0^/sec, a strain rate value used in current day state-of-the- 
art MD simulations, and 10^/sec. To the best of our knowledge, the latter is 
a strain rate several orders of magnitude slower than any reported calculation 
for a nanowire, and is a value that can actually be achieved in laboratory 
experiments on nanowires. 



3 




compressive strain e 



Figure 4: The stress versus nominal strain curve of Au nanowire under (001) compression 
(with inverted sign of stress). The data represents stress relative to surface stress at 
strain as explained in text. The stress-strain plots are shown for two different strain rates. 
The green line denotes our calculations for a strain rate of 5xl0''/sec, which is a commonly 
used strain rate in current-day MD. The red line shows our calculations for a strain rate 
of lOVsec. 

For the strain rate of 5xl0'^/sec, it was sufficient to perform ordinary 
NVT MD simulations using a time step of 2xl0~^^ sec and a Langevin ther- 
mostat with a coupling constant lxlO~^^/sec. For the strain rate of 10^/sec 
which can not be achieved through plain MD, we used our hybrid MC-MD 



algorithm as described in the preceeding section and in Tiwary and van de 



Walle (2011). The time-step and thermostat were same as for the strain rate 
of 5xl0^/sec. All the calculations were performed on our in-house parallel 
MD package. Starting value of Vq was picked such that it gave a rough 
of around 1 nanosecond (see table |l|)- This was a value high enough for our 
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current application. As the temperature/driving force (stress or strain) de- 
crease, one would pick a higher Vq that would accordingly lead to a larger 
ty^. During compression as work is being performed on the system, there is 
a change in the potential energy, and as such the value of Vq was updated 
every few thousand MD steps by an amount equaling the change in the mean 
potential energy over these many MD steps. A sharp drop in the potential 
energy indicates that a (partial) dislocation nucleation event has occurred in 
the system (see figure |5]) . Table [l] provides values of the various parameters 
used in the compression testing experiment. 



For both of these strain rates, the yielding is through slip and not twin- 
ning or elastic instabilities: a leading Shockley partial nucleates on a {111} 
slip plane at lower stresses than a trailing partial. This can be seen in figure 
|3|^b) where the leading partial nucleated from the surface and left behind 
a 2-layer thick HCP region. The trailing and leading partials are in agree- 
ment with what one expects by calculating relative Schmid factors: for (001) 
compression, (a/6) [112] and (a/6) [211] were found to be the leading and 
trailing partials corresponding to Schmid factors of 0.47 and 0.24 respec- 
tively (Weinberger et al. , 2012). 



Table 1: Starting value of Vq and expected value of for unstrained samples at various 
temperatures. For strained samples, the change in potential energy during the process of 
straining was calculated, and Vq was changed by this amount. For the temperatures 
between 350K to 425 K, ordinary MD was sufficient for the stress range considered in this 
work and hence the parameters below are only for temperatures till 325 K. 



T(K) 


Vo (eV) 




275 


-7275.00 


1400 


300 


-7267.75 


1700 


325 


-7260.25 


30 



Figure |4] brings out the perils of using unrealistic strain-rate MD calcu- 
lations as a tool for insight and discovery. We find that though the failure 
mechanism stays same for both the strain-rates, there is a significant differ- 
ence in the strain at which slip occurs. At the high strain-rate of 5xl0^/sec, 
the wire withstands strain of as high as almost 10% before the first partial 
dislocation is emitted (corresponding to a stress span of around 2.6 GPa). 
However at the more realistic strain rate of 10^/sec, the partial is emitted 
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around 8.6% only (corresponding to a stress span of around 2.35 GPa). For 
both the strain-rates, there is a distribution of the strain at which shp occurs 
and the values in Isl) denote mean values. 




Figure 5: Typical variation of the energy lid Vq (black line) and potential energy (red line) 
as a function of time. Vq is adjusted as per the change in mean potential energy every few 
thousand MD steps. A sharp drop in the potential energy and accordingly in Vq indicates 
that a nucleation event has occurred. 



This strong difference is in accordance with the strain-rate sensitivity in 
true nanowires (i.e. wires less than lOOnm in diameter) as predicted by Zhu 



et al. (2008 ) and observed in real experiments on small nanowires by Jennings 



et al. (2011). To understand and motivate this dependence, we look at the 



rate of nucleation of leading partial dislocation, as given by Eq. (10) below 



(Zhu et al, 2008): 



R = Nh'{e).exp{- 



(10) 



Here F is the Helmholtz free energy of activation as a function of tem- 
perature T and strain e (since our experimental set-up is a constant strain 
situation), ksT is the thermal energy, N is the number of equivalent surface 
nucleation sites and z/(e) is an athermal strain-dependent attempt frequency. 
Eq. ( 10 ) thus has two contributions: an athermal part related to the elastic 



limit (which we defined as stress for nucleation of the first dislocation) of the 
surface at which nucleation would occur spontaneously without any thermal 
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contributions, and an activated part that takes into account the role of ther- 
mal fluctuations in causing nucleation to happen even below the athermal 
strain (which is the minimum strain at which nucleation would occur at ab- 
solute zero temperature). 



2.4- Activation Parameters 
2.4- i- Activation Volume 

The activation volume Q is deflned as the derivative of the activation free 
energy with respect to stress, i.e. il{a,T) = — (|^)^. As reported through 
experimental measurements as well as TST based calculations, the activation 
volume for surface nucleation remains in a characteristic range of few 
(where b is the burgers vector). In comparison, for a typical bulk dislocation 
source the activation volume is upwards of 1006'^ and can be as high as 10006^ 



(Jennings et al. , 2011). The activation volume in turn determines whether a 
process will be strain-rate sensitive or insensitive. This can be reasoned as 
follows. Assuming a simple case where the activation energy depends linearly 



on stress (see Zhu et al. (2008) for detailed derivation), one can show that 
the most probable estimate of the nucleation stress is given by 

O" — CTathermal JT^^Q^ ^ ' 

where E is the Young's modulus and (Jathermai is the athermal nucleation 
stress causing instantaneous dislocation nucleation. As can be seen in Eq. 



(11), a high activation volume (as in the case of bulk dislocation source) 
masks out the effect of strain-rate. As the activation volume decreases to- 
wards values relevant for surface nucleation, the effect of strain rate should 
become very signiflcant. Figure |4] provides the flrst direct MD based evidence 
of this. 



2.4-2. Activation Free Energy 

We now report detailed calculations of the activation terms in Eq. (10). 
To do so, we performed the compression testing at a strain-rate of 10^/sec 
at 7 different temperatures from 275 K to 425 K at intervals of 25 K. These 
compression tests were stopped at various values of the nominal strain from 
8% to 8.5% (the athermal strain at Kelvin for nanowire of these dimensions 
was found to be around 13.5%). The wisdom behind choosing this particular 
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Figure 6: Strain and temperature dependence of the dislocation nucleation rate i?(e,T), 
converted here to R{a, T) by assuming a Hnear dependence of stress on strain. Each data 
point was calculated by averaging over 16 samples. The size of the individual markers 
corresponds to the 90% confidence interval in the measurement. 



range of strain will be clear soon when we provide estimates of the nucle- 
ation rate. For each of these strains, the wire is still in the elastic regime. 
As described in the introduction, we are interested in collecting statistics of 
the waiting time before nucleation of the first dislocation as the nanowire is 
held at this strain (Zuo et al. , 2005 Zuo and Ngan, 2006). The structures 
from the compression tests stopped at varying strains served as samples for 
our waiting time statistics tests. 



For each of these structures (corresponding to a combination of imposed 
strain and temperature), 16 different runs were carried out where the nanowire 
was held at a particular strain and temperature. Each of the runs was car- 
ried out until nucleation of the first dislocation, marked by appearance of a 
2-layer thick HCP region, as well as sudden dips in the potential energy and 
the stress (see figure [s]). The average rate of nucleation was then calculated 
as 

i?(e,T) = -^ (12) 

'average 
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where Taverage IS the time for nucleation averaged over the 16 samples. 

Figure [6] provides the value for nucleation rate for various strain and 
temperature combinations. In this figure we have converted strain to stress 
by assuming a roughly constant Young's Modulus of around 26 GPa as we 
obtained in figure |4] for (001) compression of Gold nanowire, in order to fa- 
cilitate comparison with published literature, even though in figure |4] we see 
a softening in the tangent modulus as the strain increases. This plot clari- 
fies our choice of imposed strains - with 8% strain (or 2.08 GPa stress), the 
nucleation rate is already slower than one every few milliseconds at 275K. 
The other end of 8.5% was picked because as illustrated in figure |4| the wire 
slips at high temperatures around 8.6%. The lengthiest of these calculations 
took around a few CPU days. With a slightly more aggressive choice of Vq, 
it should be possible to reach the one per second or still slower regime. 



For each strain e, we picked a sufficiently high value of reference tempera- 
ture To such that the rate R(e, Tq) did not any longer depend on the choice of 
temperature. We can then make the approximation that F{e, T) ^ F{e, Tq), 



and express Eq.(lO) as Eq.(13) below, to factor out the athermal frequency 



term. Given this approximation, the entropy of activation we calculate sub- 
sequently is effectively measured relative to the high temperature limit (since 
the activation entropy at high temperature could still be nonzero). 



F{e,T) ^ -knTln{§^) 



(13) 



From Eq.(13) we directly calculate the Activation Free Energy F(e,T) 
(figure 



Our values are in the rough benchmark of the 0.3 eV value 
found in Copper nano-indentation experiments where one expects homoge- 



neous nucleation to be the mechanism at work (Schuh et al. , 2005). Figure 7 



also provides the only published values of activation free energy at Kelvin 



temperature for Au nanowires (Weinberger et al. , 2012). Weinberger et al. 



(2012 )'s calculations are for a 5nm diameter nanowire (thus 4 times as many 
atoms as in our nanowire) using a chain of states methodology at Kelvin. 
A comparison between our and their free energies is thus not really justified 
due to differing system sizes - this can be understood by looking at Eq. (11 ). 



In a larger sample the number of nucleation sites (surface atoms) N is higher. 
As such the nucleation stress goes down, increasing the probability of nucle- 
ation for the same driving force (stress/strain and temperature). This in 
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turn leads to a lower free energy of activation. Some qualitative differences 
between these two calculations may will also arise because of difference in in- 



teratomic potentials. Even though a direct comparison between Weinberger 



et al. (2012)'s and our results is not justified due to these reasons, we still 



provide their results in figure [7j) since viewed together our results give a full 
picture of how the activation free energy varies with stress, temperature and 
specimen size. 




Figure 7: Activation free energy for surface nucleation of dislocations as a function of 
stress and temperature. Comparisons are made with results from literature for Kelvin 
activation energy for a wire 8 times bigger (thus with many more possible nucleation 
sites). The size of the individual markers corresponds to the 90% confidence interval in 
the measurement. 

Our calculations demonstrate how strongly and rather non-linearly the 
free energy of activation depends on the temperature. It has been a com- 
mon practice, mostly arising from lack of methods capable of providing high 
temperature activation energy barriers calculations, to assume the same tem- 
perature dependence for activation free energy across temperatures. Many 
workers have found direct and indirect evidence suggesting this is incorrect 



for studying deformation in materials. For example, Warner and Curtin 



(2009) found that in Al, a temperature dependent activation barrier can lead 
to a transition from twinning to full dislocation emission back to twinning 
with increasing temperature. We believe our algorithm should now provide 
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researchers with a tool to calculate such barriers at various temperatures 
under realistic loading rates for the first time. 



2.4-3. Activation Entropy 

From the variation of activation free energy with stress and temperature 
in figure [7| we calculated the stress dependent activation entropy. To obtain 
this quantity, we do a linear fit between the activation energy and the tem- 
perature at each stress value. The entropy is then the slope of this linear fit 
(with a negative sign), reported in figurejsj Such a calculation has rarely been 
performed for dislocation nucleation or for other problems - the two instances 
of such calculations we could find were |Hara and Li (2010)'s recent adaptive 
strain-boost hyperdynamics (ASBHD) where the authors calculate the stress 



dependent entropy for corner nucleation in Copper nanowires, and Ryu et al. 



(2011 )'s Umbrella Sampling based calculations (in Copper as well). We find 
that the entropy decreases as the driving force (stress or strain) increases, 
and is typically in the range 20-30/cb- This is roughly in the benchmark 
of values reported through previous simulations. We avoid making detailed 
comparisons here between our values and these previously published values, 
given that we differ in elements (gold versus copper), geometries (circular 
versus square with sharp cross-sections), size and ensemble (constant stress 
versus constant strain, 



see 



Ryu et al. (2011)). 
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Figure 8: Activation entropy as a function of stress for surface nucleation of dislocations. 
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3. Discussion 



3.1. Key results 

Our compressive tests on Au nanopillars show that the mechanism of de- 
formation stays same for strain rates from to lO'^/s, i.e., slip through 
nucleation of a leading Shockley partial dislocation on a {111} slip plane. 
However, the elastic limit (defined as stress for nucleation of the first disloca- 
tion) changes significantly as the strain-rate is changed. This is in qualitative 
agreement with the strain-rate dependence in nanopillars as observed in ex- 



periments (Jennings et al. , 2011 ) and also predicted by the phenomenological 



models of Zhu et al. (2008). We believe ours is the first fully atomistic calcu- 



lation reported across 6 orders of magnitudes in the strain-rate and reaching 
a realistic strain rate of 10^/s. 



The origin of activation entropy in dislocation nucleation, and related 
to it the rapid drop of activation free energy with temperature, can be at- 
tributed to thermal expansion in the material (Ryu et al. , 2011). As the 



temperature increases the expansion causes atoms to move away from each 
other making it easier for planes to shear and thus reducing the free energy 
barrier for nucleation. Our high activation entropy values show that not 
considering the temperature dependence of the activation energy can lead 
to nucleation rate being erroneous by as much as 8-12 orders of magnitude. 



This has been emphasized in the very recent literaure by Ryu et al. (2011). 



We also found that the entropy decreases as the driving force for nucleation 
(stress or strain) increases, leading to a significant and non-linear dependence 
of the activation free energy on driving force and temperature. Our method 
now provides an easy-to-implement way to calculate this dependence under 
realistic driving forces for any general sample geometry. 



3.2. Comparison of our algorithm with other algorithms for extended time- 
scales 

It is also instructive to compare our algorithm with other realistic time- 
scale algorithms for studying dislocation nucleation. Two algorithms that 
have been recently used to model temperature-dependence of the disloca- 
tion nucleation process include the Umbrella Sampling method (proposed by 
Ryu et al. ( |2011 )) and the Adaptive Strain Boost Hyperdynamics (ASBHD) 
method (proposed by Hara and Li (2010)). Both are excellent methods, but 
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have their respective hmitations. The former needs to define an order param- 
eter (or a reaction co-ordinate) for the system, in terms of which the biasing 
potential is then imposed on a certain sub-set of atoms. Picking such an or- 
der parameter might be a non-trivial task in complicated sample geometries 



and can have its pitfalls (Dickson et al. , 2009). Also, while this method is 



excellent for computing free energy barriers - it does not perform any actual 
dynamics. ASBHD does away with the need to pick up an order parameter. 
However the speed-up in ASBHD relative to ordinary MD becomes signifi- 
cant (i.e. 4-5 orders of magnitude or more) only as the temperature of the 
system falls below lOOK or so. The method is a local boost scheme, thus 
specific atoms are lifted out of the low-lying energy basins, making them 
preferential sites for nucleation to happen. Such a local boost scheme works 
well for specific geometries (such as a square nano-rod with sharp corners), 
but might not be very well suited for studying more homogeneous nucleation 
as we are considering in the current paper. Since our method does not require 
the specification of the degrees of freedom of interest (which is crucial when 
the mechanisms are complex and involve the movement of many atoms), it 
is well suited for studying homogeneous nucleation. 



By providing a time-scale correction independent to the main simula- 
tion, our algorithm also provides the ability to implement time-dependent 
boundary conditions (relevant to a tensile test for e.g.) and in general time- 
dependent forces. To do this, we can launch a set of adiabatic switching jobs 
prior to the main simulation, that give a good starting value for the quan- 
tity ti^VLS for the boosted time-scale. This is to be contrasted with 



ASBHD (Hara and Li 2010) and Hyperdynamics methods in general (Voter 



1997) where time-scale estimates remain noisy and non-converged for long 



simulation times, especially as one tries to increase the speed-up relative to 
ordinary MD (Miron and Fichthorn 2003 Tiwary and van de Walle 2011). 



4. Conclusions 

In this paper we have derived and demonstrated a hybrid MC-MD al- 
gorithm that can be used to achieve realistic time-scales in fully atomistic 
simulations of materials while still predicting correct deformation physics. 
The algorithm is especially designed to be suited for massive parallelization. 
By using this algorithm, we obtained compression testing stress-strain plots 
at strain rates several orders of magnitude lower than ever previously re- 
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ported for MD simulations. We showed that high strain-rates in simulations, 
which have been common due to lack of methods capable of implementing 
low strain-rates, can lead to a significant error in the elastic limit (defined 
earlier) of the material. We also derived the full stress and temperature de- 
pendence of the activation free energy for surface nucleation of dislocations 
in Gold nanowires. The algorithm was described in sufficient detail to be 
useful to the mechanics community for different applications. 
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